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Spectral correlations in unitary invariant, non-Gaussian ensembles of large random matrices pos- 
sessing an eigenvalue gap are studied within the framework of the orthogonal polynomial technique. 
Both local and global characteristics of spectra are directly reconstructed from the recurrence equa- 
tion for orthogonal polynomials associated with a given random matrix ensemble. It is established 
that an eigenvalue gap does not affect the local eigenvalue correlations which follow the universal 
sine and the universal multicritical laws in the bulk and soft-edge scaling limits, respectively. By 
contrast, global smoothed eigenvalue correlations do reflect the presence of a gap, and are shown to 
satisfy a new universal law exhibiting a sharp dependence on the odd/even dimension of random 
matrices whose spectra are bounded. In the case of unbounded spectrum, the corresponding uni- 
versal 'density-density' correlator is conjectured to be generic for chaotic systems with a forbidden 
gap and broken time reversal symmetry. 
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I. INTRODUCTION 

Ensembles of large random matrices H gener- 
ated by the joint distribution function P [H] oc 
exp{— /3Tr V [H]}, with (3 being a symmetry parameter 
as explained below, may display phase transitions under 
non-monotonic deformation of the confinement potential 
V [H] . Different phases are characterized by topologically 
different arrangements of eigenvalues in random matrix 
spectra that may have multiple-band structure. Ran- 
dom matrices, whose spectra undergo phase transitions, 
appear in quantizing two-dimensional gravity in 
the context of quantum chromodynamics QJ^], as well 
as in some models of particles interacting in high dimen- 
sions . Transition regimes realized in invariant random 
matrix ensembles have implications for a certain class 
of Calogero-Sutherland-Moser models |7j] . These matrix 
models may also be applicable to chaotic systems having 
a forbidden gap in the energy spectrum. 

In the eigenvalue representation, the invariant random 
matrix model is defined by the joint probability distribu- 
tion function @1 



JV 



JV 



P({e})=Z N 1 n ki-e/n^H 37 ^ w 



i>j=l 



fc=l 



of N eigenvalues {e} = {ex, ejv} of an N x N random 
matrix H. The symmetry parameter f3 coincides with a 
number of independent elements in oft-diagonal entries 
of a random matrix H. For real symmetric matrices, 
(3 = 1 (orthogonal symmetry), [3 — 2 for Hermitian ma- 



trices (unitary symmetry), and (3 = 4 for self-dual Her- 
mitian matrices (symplectic symmetry). It is convenient 
to parametrize the confinement potential V (e) entering 
Eq. (pi) by a set of coupling constants {d} — {d\, d p }, 



V{e) = j2^/\d p >Q, 

k=l 



(2) 



so that we may consider the phase transitions as occur- 
ring in {d}-space. Because the confinement potential is 
an even function, the associated random matrix model 
possesses so-called Z2-symmetry. 

Variations of the coupling constants affect the Dyson 
density vd, that can be found by minimizing the free en- 
ergy Tn = ~ log-Ziv, Eq. (|l|), subject to a normalization 
constraint J vr> (e) de = N, 



dV 

de 



-V 



dC-^r = 0, 



(3) 



where V indicates a principal value of the integral. When 
all dk are positive, so that confinement potential is mono- 
tonic, the spectral density v D has a single-band support, 
Ab = 1. Non-monotonic deformation of the confinement 
potential can be carried out by changing the signs of 
some of dk (k ^ p). Such a continuous variation of cou- 
pling constants may lead, under certain conditions, to a 
discontinuous change of the topological structure of spec- 
tral density vjj, when the eigenvalues {e} are arranged 
in A/b > 1 'allowed' bands separated by 'forbidden' gaps. 

The phase structure of Hermitian ({3 — 2) one-matrix 
model Eq. (Q) has been studied in a number of works 
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12j, where the simplest examples of non-monotonic 
quartic and sextic confinement potentials have been ex- 
amined, ft has been found that there are domains in the 
phase space of coupling constants where only a particu- 
lar solution for Vn exists, and it has a fixed number At, 
of allowed bands. However, in some regions of the phase 
space, one can have more than one kind of solution of the 
saddle-point equation Eq. (|J). In this situation, solutions 
with different number of bands A/^ , ... are present 

simultaneously. When such an overlap appears, one of 
the solutions, say A/^ , has the lowest free energy 
and this solution is dominant, while the others are sub- 
dominant. Moreover, numerical calculations |l2| showed 
that some special regimes exist in which the bulk spectral 
density obtained as a solution to the saddle-point equa- 
tion Eq. (||) differs significantly from the genuine level 
density computed numerically within the framework of 
the orthogonal polynomial technique. It was then ar- 
gued that such a genuine density of levels cannot be in- 
terpreted as a multi-band solution with an integer num- 
ber of bands. A full understanding of this phenomenon 
is still absent. 

Recently, interest was renewed in multi-band regimes 
in invariant random matrix ensembles. An analysis based 
on a loop equation technique 13 lj] showed that finger- 
prints of phase transitions appear not only in the Dyson 
density but also in the (universal) wide-range eigenvalue 
correlators, which in the multi-band phases differ from 
those known in the single-band phase 1 15 - 17 1 . A renor- 
malization group approach developed in Ref. |Tg| ] sup- 
ported the results found in Refs. |L3 14 for the partic- 
ular case of two allowed bands, referring a new type of 
universal wide-range eigenlevel correlators to an addi- 
tional attractive fixed point of a renormalization group 
transformation. 

The method of loop equations p3|,^4[ , used for a treat- 
ment of non-Gaussian, unitary invariant, random ma- 
trix ensembles fallen in a multi-band phase, is only suit- 
able for computing the global characteristics of spectrum. 
Therefore, an appropriate approach is needed capable of 
analyzing local characteristics of spectrum (manifested 
on the scale of a few eigenlevels) . A possibility to probe 
the local properties of eigenspectrum is offered by the 
method of orthogonal polynomials. A step in this direc- 
tion was taken in a recent paper |l9| l , where an ansatz was 
proposed for large-A^ asymptotes of orthogonal polyno- 
mials associated with a random matrix ensemble having 
two allowed bands in its spectrum. Because the asymp- 
totic formula proposed there is of the Plancherel-Rotach 
type p0[ , it is only applicable for studying eigenvalue 
correlations in the spectrum bulk and cannot be used 
for studying local correlations in an arbitrary spectrum 
range (for example, near the edges of two-band eigen- 
value support). 

The aim of the present paper is to develop a new ap- 



proach (within an orthogonal polynomial scheme) allow- 
ing a unified treatment of eigenlevel correlations in the 
unitary invariant U (N) matrix model {(3 = 2) with a for- 
bidden gap. This is a further extension of the Shohat 
method ]2l],^2| that has been used previously by the au- 
thors to study U (N) invariant ensembles of large random 
matrices in the single-band phase ]23|,Q. In particular, 
we are able to study both the fine structure of local char- 
acteristics of the spectrum in different scaling limits and 
smoothed global spectral correlations. Our treatment is 
based on the direct reconstruction of spectral correla- 
tions from the recurrence equation for the corresponding 
orthogonal polynomials. 



II. GENERAL RELATIONS 

In this section we briefly review the orthogonal polyno- 
mial technique (|] . The n-point correlation function that 
describes the probability density to find n levels around 
each of the points Si, e n when the positions of the re- 
maining levels are unobserved is defined by the formula 



Nl 



(N — n)\ 



/ P({e}) de k . (4) 

J — OO j. in 



This correlation function can explicitly be expressed in 
terms of the two-point kernel (e, e') as follows 



R n (ex, ...,£„) = det \\K N {e h e 



J ) \\i,j=l...n 



(5) 



Here, 



v , t\ <Pn (gQ fN-i (e) - (e) <Pn-i (gQ 

K N {£,£) = C N ; , (6) 

e' — e 

and the 'eigenfunctions' 

<p n (e) = P n (e)eiq>{-V(e)} (7) 

are determined by the set of polynomials orthogonal with 
respect to the measure dfj, (e) = exp {—2V (e)}de, 



d\i (e) P n (e) P m (e) = S nm , 



and obeying the recurrence equation 

eP„_i (e) = c n P n (e) + c n ^P n - 2 (e) . (9) 

The recurrence coefficients c„ entering Eqs. (|^) and (|^) 
are uniquely determined by the measure dfi. Equations 
(||) and (^) demonstrate that the problem of eigenvalue 
correlations is reduced to that of finding asymptotes for 
the eigenfunctions (fiN- 
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III. MAPPING RECURRENCE EQUATION 
ONTO DIFFERENTIAL EQUATION 



A. Single— band phase 



To map a recurrence Eq. (||) onto a second-order dif- 
ferential equation for cigcnfunctions ip n , we note that the 
first derivative dP n /de can be represented as pl|,p2| 



dP 

^=A n (e)P n _ 1 -B n (e)P n 
de 



where 

A n (e) = 2c„ / dfj, (t) 



(10) 



(11) 



B n (e) = 2c n J dn (t) 



V (t)-V> (e) 



P n (t) P n _ x (t) . 



(12) 

Then, by using Eqs. @ and @, one obtains after some 
algebra that the fictitious wave function tp n given by Eq. 
(0) satisfies the following differential equation: 

d 2 Lp n (e) dip n (e) 

- .F n ye) — h Q n (e) ip n (e) = 0. (13) 



de 2 



de 



Here. 



and 



Gn (e) 



T n (e) 

dB n ^ _c, 
de 

—B n 
d 2 V 



Cn-l 



de 2 



1 dA n 

A n de 



-A n A n - 

2— 

de 

dV_ x 2 

de 



(14) 



1 dA n 
A n de 

1 dA n dV 



A n de de ' ^ 

Equation ( |l3| ) is valid for arbitrary n. We note that de- 
spite the generality of the differential equation obtained, 
its practical use is quite restricted since the functions 
T n (A) and Q n (A) entering Eqs. Jl3l ) can be calculated 
explicitly only for rather simple measures d/i. Neverthe- 
less, an asymptotic analysis of this equation is available 
in the limit n — N ^> 1, that is of great interest in ran- 
dom matrix theory. 



The single-band phase corresponds to monotonic con- 
finement potentials or to those having light local extrema. 
Corresponding asymptotic analysis has been carried out 
by the authors in Refs. [ p3]]24j| . For further compari- 
son with a two-band-phase solution, we give a differen- 
tial equation for tpjj (e) obtained in the leading order in 
N 3> 1 [upper index indicates that the single-band phase 
is considered]: 



d 2 ^ (e) 
de 2 




d<P% (e) 
de 



,(!) 



00 



(16) 



It is remarkable that Eq. ( |l6| ) does not contain the con- 
finement potential explicitly, but only involves the Dyson 
density 



tP 



'Pa 



tdt dV 1 



e 2 /V 2 



N 



f 2 



dt 



l-t 2 /V 2 N 



(17) 



corresponding to the single-band phase and analytically 
continued on the entire real axis; T>n is the soft edge 
of the spectrum, being the positive root of the integral 
equation 



dV 



tdt 



dt ^/Vj~^¥ 



~2~' 



(18) 



It has been shown that for non-singular confinement po- 
tential, solutions of Eq. ( |l6|) lead to the universal sine 
kernel in the bulk scaling limit, and to the so-called G~ 
multicritical correlations in the soft-edge scaling limit 
pjj . An additional logarithmic singularity of confine- 
ment potential introduces additional terms into Eq. (^) , 
giving rise to the universal Bessel correlations in the ori- 
gin scaling limit p5|]23j|. F or further progress in the field, 
see very recent paper |26[ . 



B. Two— band phase 



Let us consider the situation when the confinement potential has two deep wells leading to the Dyson density 
supported on two disjoint intervals located symmetrically about the origin, T>~^ < \e\ < Pi. In this situation, the 
recurrence coefficients c„ entering Eq. (^J) are known to be double-valued functions of the number n [p]JTo| . This 
means that for n = JV ^> 1, one must distinguish between coefficients CM±2q ~ cn and coefficients c/v-i±2g ~ C/v-ij 
belonging to two different smooth (in index) sub-sequences; here, integer q ~ O [N J . Bearing this in mind, the 
large-TV version of recurrence equation Eq. (p|) can be rewritten as 
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[e 2 - (c 2 N + c 2 v _ 1 )] P N (s) = c N c N -! [P N -i (e) + P N+1 (e)] , 
whence we get the following asymptotic identities: 



e 2 ^( £ ) = (4 + 4_ 1 ) A E 



k=0 



C N + C 7V-1 



k k 

E 

3=0 



PN+4,j-2k (e) 



(19) 



(20) 



and 



, 2A+1 P.v(e) = (4 + 4_ 1 ) A E 



fe=0 



CATCtv-1 



-N ~ ^N-l 



k k 

E 

3=0 



[cN-lPN+4j-2k+l (e) + C7VpV+4j-2fc-l (ff)] (21) 



with integer A > 0. 

Expansions Eqs. ( p0| ) and ( pl| ) make it possible to compute the required functions J-jv and (?jv entering the 
differential equation Eq. ( |l3| ) for fictitious wave functions in the limit N 3> 1. Substituting the explicit form of the 
confinement potential set by Eq. (0) into Eqs. (Kw and (O), we obtain 



(22) 



p 2fc-l „ 

A N {e) = 2c N Y^d k E e* -1 / dfi(t) P% (i)* 2 ^" 1 , 



and 



fe=l A=l 



p 2fe-l 



Bjv (e) = 2c w E dfe E £A_1 / ^ (*) Pw (*) P ^-! W 1 



2&-A-1 



(23) 



k=l A=l 



respectively. Both integrals above can be calculated using expansions Eqs. (0), (pl|), and exploiting the orthogonality 
expressed by Eq. (||). Detailed calculations, given in Appendices A and B, lead to the following results: 



A N (e)^^{v+-(-l) N V N )v£ 



« dV t 2 



dt 



dt t 2 - e 2 



(V+Y-t 2 t 2 -(V- N )' 



(24) 



B N (e) = -eV 



v - dV 



i 2 -{-ifv- N v N 



dt 



dt 



(v + N f-t 2 t 2 -{v- N y 



dV 



t 2 - e 2 de 



(25) 



Having obtained the explicit expressions for functions An and Bn, it is easy to verify that coefficients T n (e) and 
Q n (e) entering the differential equation Eq. (^|) for the fictitious wave function ipn (e) may be expressed in terms of 
the Dyson density in the two-cut phase supported on two disconnected intervals e G (— "D%, ^^n) U (^n^n) 



V 



N) 



V 



dt 



dV/dt 



t 2 -e 2 



when N > 1. Namely, Eqs. §L§, (|15|), (g§ and (]25j) yield 



(e) = -j- lo. 



7r| e |^ n) (e) 



- to 



(26) 



V 



N! 



(27) 



Sat (e) 



,(«) 



,(«) 



Nl 



- e z 



e 2 -{v N y 



s 2 + (-l) N V- N V + N 



(28) 
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In the large-iV limit, the second term in Eq. ( |28| ) can be neglected provided e belongs to the one of allowed bands, 
so that ipffi (s) satisfies the following asymptotic differential equation in the two-cut phase: 



de 2 



( 



de 



log 



n\e\u^ ( £ ) 



V 



N) 



V 



df [ N ] (g) 
de 



(II) / si (II) / \ n 

tti/), (e) <p N (e) = 0. 



(29) 



We recall that T> N and T> N are the end points of the eigenvalue support that obey the two integral equations 



v - dV 



t 2 dt 



dt 



t 2 - (V 



N) 



nN 
~2~' 



(30) 



dt 



dt 



= 0, 



- V 



(31) 



obtained in Appendix C. One can verify that as T> N tends to zero, we recover equation Eq. ( |16[ ) valid in the single-band 
regime. 



IV. LOCAL EIGENVALUE CORRELATIONS 



Eigenvalue correlations in the spectra of two-band random matrices are completely determined by the Dyson density 
of states entering the effective Schrodinger equation Eq. (p9|). 

(i) In the spectrum bulk, the Dyson density is a well-behaved function that can be taken approximately as a 
constant on the scale of a few eigenlevels. Then, in the vicinity of some £q that is chosen to be far enough from the 



spectrum end points ±2?^, Eq. (29) takes the form 



d 2 ^' (e) 
de 2 



+ [7r/A( £0 )r^ I) (e)=0, 



(32) 



with A (eq) = (eo) being the mean level spacing in the vicinity of eq. Clearly, the universal sine law for the 

two-point kernel, Eq. (^|), follows immediately. 

(ii) Eigenvalue correlations near the end points of an eigenvalue support are determined by the Dyson density as 
well. Noting that in the vicinity of |e| = 2?^, the Dyson density can be represented in the form p7|j2~^ ]. 



(ii) 



± 1 



V 



(33) 



where TZn (il) ^ and m is the order of multicriticality, we readily recover the universal multicritical correlations 
previously found |2~f| in the soft -edge scaling limit for U (N) invariant matrix model in the single-band phase. 



V. SMOOTHED CONNECTED 'DENSITY-DENSITY' CORRELATOR 



Let us turn to the study of the connected 'density-density' correlator that is expressed in terms of the two-point 
kernel, Eq. (||), as follows 

(Su N (e) 8v N (e')) u = - ° N {<& (e) <^ 2 V _ 1 (e') + <p 2 N (e') ^ 2 N _ X (e) - 2<p N (e) <p N ^ (e) tp N (e r ) (e')} (34) 
(s - e') 

where e / e', and the upper index (II) in (p n has been omitted for brevity. We still deal with the two-band phase. 
Expression Eq. d34) contains rapid oscillations on the scale of the mean level spacing. These oscillations are due to 
presence in Eq. (p3) of oscillating wave functions ipN and (fN-x- 
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To average over the rapid oscillations, we integrate, over the entire real axis, rapidly varying wave functions in Eq. 
multiplied by an arbitrary, smooth, slowly varying function. To illustrate the idea, consider the integral 



If- ' -l- 2 



de<p J N (s)f(e), (35) 
where f (e) is arbitrary slowly varying function that should be chosen to be even due to the evenness of (p 2 N (e). Setting 

oo 

f(e) = ^f Q e 2Q , (36) 



we immediately obtain with the help of Eqs. (Al) and (A6) that 



^ 2 r v * 

If = V f Q A 2Q = - / 

IT IT>- 

a=0 N 



ef (e) de 



(KY 



V 



N) 



(37) 



Bearing in mind that both f (e) and <p 2 N (e) are even functions, the last integral can be transformed as follows 



If = 



+ 00 



de<p% (e) f (e) = - 



\e\ f (e) de 



V 



N J 



(38) 



whence we conclude that in the large- N limit, 

2 / x 1 lei 

<Pn ( £ ) = ; 



(KY 



V 



N) 



Q{V+-\e\)Q{\e\-V 



Nl 



(39) 



The same procedure should be carried out with expression tp^ (e) <Pn—i ( £ ) m Eq. (p4|). Since this construction is 
an odd function of e, we have to consider the integral 



I g = dsip N (e)<^Ar_i (e)g(e) , 



with 



Q=0 

being a smooth odd function. It is easy to see with the help of Eqs. (Bl), (B5) and (C7) that 



(40) 



(41) 



Ig — X/ S" r 2 Q +l — 



a=0 



v+-{-ifv- 



vt g( £ ) e*-{-l) N V- N V N 



ds 



V 



N) 



e 2 - V 



N / 



(42) 



Exploiting the oddness of g (e) and <pn (s) <PN—i (s), we write Eq. ( f42| ) in the form 

r+oo ^ 



deip N (e) ip N ^i (e)g(e) 



N 



e 2 -(-!)" V N V N 



sgn (e) de 



JV I I JV 



V 



N) 



e 2 - V 



N) 



(43) 



Equation (Ej3) leads us to the conclusion that in the large- A limit, 



(f N (e) iptt-x (e) 



sgn (e) 



e*- { -lfv- N V + N 



N 



V 



Nl 



e 2 - (V 



N) 



®{V+-\e\)Q{\e\-V 



N) 



(44) 



Combining Eqs. (|34|), (p9|), (|4J) and (C7), we finally arrive at the following formula for smoothed 'density-density' 
correlator, 



(6u N (e) Su N (e'))„ = (V+ \e\) 9 (|e| - V N ) 9 (V+ - \e'\) 9 (|e'| - V N ) 



1 



(e-e'Y 



ze' - (V 



N) 



(V 



Nl 



EE 



[V 



N) 



(V 



Nl 



v 



N) 



(45) 



V 



Nl 



+ (-1) 



N 



IK) 



- v 



Nl 



The same formula can be obtained by WKB by solving Eq. (|29|), using definition Eq. (M) followed by averaging over 
rapid oscillations. It can be verified that for N even, this result coincides with Eq. (6.6) of Ref. |Q where it was 
obtained by a completely different method, and for the case of odd N being omitted. 

It is seen from Eq. ( f45| ) that smoothed 'density-density' correlator in the two-band phase is a new universal 
function in random matrix theory. It is universal in the sense that the information of the distribution Eq. ([!]) is 
encoded into the 'density-density' correlator only through the end points T>^ of the eigenvalue support. A striking 
feature of the new universal function Eq. ( ft5| ) is its sharp dependence on the oddness/evenness of the dimension N 
of the random matrices whose spectra are bounded. The origin of this unusual large-TV behavior will be discussed in 
the next Section. 



Finally, let us speculate about the universal correlator Eq. 



in the limit of unbounded spectrum, T>J, — > oo, with 



a gap. Inasmuch as it describes correlations between the eigenlevels which are repelled from each other in accordance 
with the logarithmic law, that is known to be realized (2^J|^] in the weakly disordered systems on the energy scale 
\e — e'\ <C E c (E c is the Thouless energy), we may conjecture that the corresponding limiting universal expression 



lim {8v N (e) 8v N (e')) n = - n , 
v+-*+oo 2ir 2 (e 



sgn (ee') 



ee 



e'Y Vie 2 ~ A2 ] [e' 2 - A 'l 



9 (|e| - A) 9 (|e'| - A) 



(46) 



reflects the universal properties of real chaotic systems with a forbidden gap A = T>Jj and broken time reversal 
symmetry, provided \e — e'\ <C E c , In two limiting situations (i) of gapless spectrum, A = 0, and (ii) far from the gap, 
e| , |e'| ^> A, the correlator Eq. ( (46| ) coincides with that known in the random matrix theory of gapless ensembles 
Jui| . fj"7| and derived in Ref. ps| ] within the framework of diagrammatic technique for spectrum of electron in a random 
impurity potential. 



VI. CONCLUDING REMARKS 



In this study we developed a unified formalism allowing the computation of both global and local spectral char- 
acteristics of U (N) invariant ensembles of large random matrices possessing Z2-symmetry, and deformed in such a 
way that their spectra contain a forbidden gap. We proved that in the pure two-band phase, the local eigenvalue 
correlations are insensitive to this deformation both in the bulk and soft-edge scaling limits. In contrast, global 
smoothed eigenvalue correlations in the two-band phase differ drastically from those in the single-band phase, and 
generically satisfy a new universal law, Eq. ([45|), which is unusually sensitive to the oddness/evenness of the random 
matrix dimension if the spectrum support is bounded. On the formal level, this sensitivity is a direct consequence of 
the 'period-two' behavior fj|,|l0| of the recurrence coefficients c n [see Eq. (||)] that is characteristic of two-band phase 
of reduced Hermitian matrix model. To see this, consider the simplest connected correlator (Tr HTrH) that can be 
exactly represented in terms of recurrence coefficients for any n, 

(TrHTrH) c = 4. (47) 
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Since in the two-band phase c n is a double-valued function of index n, alternating between two different functions 
as n goes from odd to even, the large- N limit of the correlator (TrHTrH) c strongly depends on whether infinity is 
approached through odd or even N. Then, an implementation of a double-valued behavior of c n into the higher order 
correlators of the form (TrH' 1 ' TrH ! ) contributing to the connected 'density-density' correlator gives rise to the new 
universal expression Eq. (J4q). 



Let us, however, point out that no such sensitivity has been detected in a number of previous studies jL3 14 
exploiting a loop-equation technique. This is due to the following reasons. In the method of loop equations, used for 
a treatment of non-Gaussian random matrix ensembles fallen in a multi-band phase, one is forced to keep the most 
general (non-symmetric) confinement potential V (e) — J2kLi dk£ k /k until very end of the calculations, leading to a 
necessity to take the thermodynamic limit N — > oo prior to any others. Therefore, Z2-symmetry in this calculational 
scheme can only be implemented by restoring Z2-symmetry at the final stage of the calculations, setting all the extra 



coupling constants d,2k+i to zero. Doing so, one arrives at the results reported in Refs. [[13 14 . 

^From this point of view, the formalism developed in this paper corresponds to the opposite sequence of thermo- 
dynamic and Z2-symmetry limits, since we have considered the random matrix model that possesses Z2-symmetry 
from the beginning. Qualitatively different large-iV behavior of the smoothed connected 'density-density' correlator, 
Eq. ([45|), and of the smoothed connected two-point Green's function given by Eq. (15) of Ref. H] provides a direct 
evidence that the order of thermodynamic and Z2-symmetry limits is indeed important when studying global spectral 
characteristics of multi-band random matrices. 
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Appendix A: Calculation of the function Am (e) 



Let us consider an integral 



A 



2a 



dn (t) P N (t) t 2 



with integer a > 0. Making use of Eq. (|20|), we rewrite A 2fT in the form 

k k 



a 2ct = (c?N+4f-iYJ2(l) 

I — n V / 



a \ I cjvCjv_i 



C N + C N-1 



E ( ] ) / dn{t)P N {t)P N+ ^ 2k (t). 



k=0 v 7 x JV -JV-i/ J=0 

Orthogonality of the P„ allows us to integrate over the measure d/i, thus simplifying Eq. (A2) 

k k 

"2j, 



a 2ct = (c 2 N + c 2 N _ 1 y^2 



k=0 



a \ I ctvCat_i 



> — ^- , E ( * ) 



where 8\, is the Kronecker symbol. Using integral representation 



S k k , = Re / — cxp {i (k - k') 6} , 
Jo 2n 

one can perform the double summation in Eq. (A3): 

A 2ct = / — (c 2 N + c^r_ 1 + 2cArcjv-i cos 6) a . 

JO ^7T 

Introducing a new integration variable t 2 = c 2 N + c 2 N _ 1 + 2cncn-i cos#, we derive an integral formula 

t 2 ° +1 dt 



A 



2 f v » 



2a 



7T JV 



(v+y-t 2 t 2 -{v- N y 



with 



V N = \c N ± cat_ 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



Now, taking into account representation Eq. (A6) for A 2(7 , and using the fact that A 2(J+ i = 0, we obtain from Eq. 



A N (e) = 2c N E dk E A 2(k-a)£ 2a ' 



(A8) 



k=l (7=1 



Summing over cr yields 



A N (e) = > d k I 

77 /x>- 



fe=i 



t 2 - (V 



t 2 -e 2 ' 



(A9) 



whence we get, with the help of Eq. 



A N (e) = 



dt 



t ( dV dV 

t— £ 



IK) -* 2 



t 2 -(v- N y 



t 2 - e 2 \ dt de 



(A10) 



Further, noting that 



dt 



t 2 - (v 



N) 



t 2 -e 2 



0, 



(All) 



and taking into account Eq. (C7), leads to the final expression given by Eq. (E 



Appendix B: Calculation of the function Bn (e) 



Let us consider an integral 



r 2ff+ i = J dfi (t) P N (t) P N - X (t) t 



2(7+1 



(Bl) 



with integer a > 0. Making use of expansion Eq. (^l|), we rewrite Eq. (Bl) in the form that allows us to perform the 
integration over the measure dfj,: 



iWi = ~ (4 + 4-0 " / d^ (t) p N ^ (t) J2 

z J fc=0 



k k y , 

a \ cnCn-i \ \— * / K 



k I \c 2 +c 2 



E 



x [CjV-l-PiV+4i-2fc+l (*) + c iV-P/V+4j-2fc-l (*)] ■ 



(B2) 



After integration, we get 
1 



r 2o -+i - - {c~ N + c N _ 



-0 ff E 



a \ I cnCn-i 



-N ' "W- 



;) fe E(j) [cN-x8% +1 + c N 8l]. (B3) 



The double summation in Eq. (B3) can be performed using the integral representation for the Kronecker symbol 
given by Eq. (A4): 



IV+i = - / — (c% + c%_ t + 2c n cn-i cosd)" [cn + CJV-1 cos#] . 

2 _ „2 „2 



Introducing a new integration variable i = c N + cj^_ 1 + 2cncn-i cos6, we get 



(B4) 



T2o-+l 



1 [< 



"""CAT J-p 



t 2 ° +1 dt 



(V+f-t 2 t 2 -{V~ N ) 



[t +c N — c JV _ 1 ] 



Then, Eqs. (ph, (Bl) and (B5) yield 



(B5) 



k-l 



Bn (e) — — dk ^ r 2 fe-2o--i£ 2<T 1 - 



fe=l er=l 

Summing over cr leads to the integral expression 

P rV+ r+2 , J2 _ „2 



^ (e) = - E d * 



It + c w Cjv— lj 



et2fc-i _ fe 2*-i 



fc=l "^N 



2 /"^ 



(^) 2 -* 2 l V 2 -(P- N ) 2 



t 2 



dt 



t 2 - v 



N) 



t 2 +c 2 N -c%_ 1 / dV^ dV 

~ V 2 vltt _ 1 ~d7 



(B6) 



(B7) 



Now, taking into account Eqs. (All) and (C6), we obtain Eq. (p5|). 
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Appendix C: Soft edges of eigenvalue support 



To find the equations determining the points T> N where the Dyson spectral density goes to zero, we start with the 
following formula from the theory of orthogonal polynomials pffl 



// -2,;, / d f i(t)^P n (t)P n - 1 {t). 



(CI) 



Let us use expansion Eq. (g^Tj) to calculate asymptotically the integral entering Eq. (CI) in the limit n — N ^> 1. It 
is easy to see that 



p p v 
N = 2c N J2dx / dfi (i) P N (t) P N ^ (t) t 2 *- 1 = 2c N d x T 2 x-i, 



\=i " x=i 
where T2A-i is given by Eq. (B5). Then, we immediately obtain the relationship 

r 

IT>Z 



N= 2 - 

7T 



dt 



t 2 - v 



dV 
' dt 



This result rewritten for n = N — 1, yields in the large- TV limit, 

r 



N= 2 - 

IT 



dt 



dV 



t 2 - V 



NJ 



dt 



Equations (C3) and (C4) yield two equations whose solutions determine the edge points U^: 



t 2 dt 



t 2 - (v 



NJ 



dV _ nN 



and 



dt 



dt 



dV 



t 2 (VnY 



dt 



= 0. 



(C2) 



(C3) 



(C4) 



(C5) 



(C6) 



Finally, we note that because P_i (e) = 0, it follows from 
Eq. ([)[) that cq = 0, and as a consequence, an even 
branch C2N always lies lower than an odd branch C2N±i, 
so that C2N < C2N±i- Then, we may conclude from Eq. 
(A7) that 



V%-(-l) N V 



cn 



N 



(C7) 
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